{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "##############################################################################\n",
    "# Script used in the paper:\n",
    "# Dune migration and volume change from airborne LiDAR, terrestrial LiDAR \n",
    "# and Structure from Motion--Multi View Stereo\n",
    "# by\n",
    "# Carlos H. Grohmann et al - 2019/2020\n",
    "# guano (at) usp (dot) br\n",
    "# Institute of Energy and Environment - University of Sao Paulo\n",
    "#\n",
    "# Please check the GitHub repo for the final reference to the paper\n",
    "##############################################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Monte Carlo analysis of random points\n",
    "\n",
    "To compare the elevation of elevation datasets, a few key points must be considered beforehand, such as spatial resolution, presence of voids and number of samples used in the analysis. A direct comparison of the raster datasets on a pixel-by-pixel basis, is not the best approach, since differences in pixel size and presence of voids could affect correlation analysis, and the large number of samples would impact descriptive statistics, goodness-of-fit and error measurements. \n",
    "\n",
    "A random sample of elevation values can overcome these issues, but raises the question of how many data points are needed to properly represent the original dataset. To answer this question, a Monte Carlo approach was devised in the following form:\n",
    "\n",
    "- the model was run 50 times;\n",
    "- the number of random points analyzed was n=50, 100, 250, 500, 1000, 2500, 5000 and 10000;\n",
    "- in each run, *n* random points were created and elevation was extracted from SRTM;\n",
    "- after 50 runs, correlation was calculated between the first and the 49 subsequent sets of *n* random points;\n",
    "- a Four Parameter Logistic Regression (4PL) was calculated for the mean value of correlation of each set of *n* random points.\n",
    "\n",
    "In order to ensure reproducibility of the analysis, the random seed used to generate points was set to the sequential number of each model run (0,1,2,3,...,49) multiplied by a constant. Results of this approach have shown that for the TLS aurvey area, 1000 random points can be used to describe the elevation of the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1 - Extract elevation for sets of randon points [50,100,250,500,1000,2000,5000,10000] and then\n",
    "# calculate stats for each. \n",
    "\n",
    "# MonteCarlo-like analysis:\n",
    "# 0 - for a series of n random points [50,100,250,500,1000,2000,5000,10000]:\n",
    "# 1 - get X sets of n random points (rand_n_01, rand_n_02, rand_n_03,...) - sorted\n",
    "# 2 - calculate correlation between first set and all others\n",
    "# 3 - put data in a table and plot the results\n",
    "# 4 - make plot off all values (X = n_points, Y = correlation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import python libraries\n",
    "import sys, os, itertools\n",
    "import numpy as np\n",
    "import scipy.stats as ss\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import subprocess\n",
    "from IPython.display import Image # can use this to display GRASS maps\n",
    "# stats\n",
    "from statsmodels.graphics.gofplots import qqplot\n",
    "from scipy.stats import linregress\n",
    "from scipy.optimize import leastsq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# helper func: delete all *random maps\n",
    "def clean_rand():\n",
    "    grass.run_command('g.remove', type='vector', pattern='*random*', flags='f')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# helper func: round to nearest 5\n",
    "def round5(x):\n",
    "    rounded = int(round(x/5.0)*5.0)\n",
    "    return rounded"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# aux func: sample DEMs elevations\n",
    "# requires dbsae connection in GRASS to be SQLITE\n",
    "def sample_dems_mc(dem, n_points, mc_run, ow_vector, ow_what, vmask):\n",
    "    ''' create random points for selected tile and sample the elevation values \n",
    "    simplified version of sample_dems tailored for MonteCarlo-like analysis\n",
    "    dem = raster map (elevation)\n",
    "    n_points = number of random points\n",
    "    mc_run = number of times a set of random pints will be created \n",
    "    ow_vector = should vector maps be overwritten?\n",
    "    ow_what = re-run v.what.rast ? \n",
    "    vmask = vector mask to restrict points\n",
    "\n",
    "    note: to keep random points really randon and yet ensure reproducibility,\n",
    "    random seed is set to the value of mc_run * 42'''\n",
    "    \n",
    "    grass.run_command('g.region', raster=dem, flags='a')\n",
    "    \n",
    "    # random points \n",
    "    vector_name = dem.split('_')[0] +'_random_' + str(n_points) + '_' + str(mc_run).zfill(2)\n",
    "    grass.run_command('v.random', output=vector_name, npoints=n_points, restrict=vmask, seed=mc_run*42, quiet=True, overwrite=ow_vector) \n",
    "    rand_col = 'rand_' + str(n_points) + '_' + str(mc_run)\n",
    "    grass.run_command('v.db.addtable', map=vector_name, columns=rand_col+' double precision', quiet=True, overwrite=ow_vector)\n",
    "        \n",
    "    # sample raster map - force if overwrite vector is true\n",
    "    if ow_vector or ow_what:\n",
    "        grass.run_command('v.what.rast', map=vector_name, raster=dem, column=rand_col, quiet=True)\n",
    "        \n",
    "    # export as ascii and read into python\n",
    "    xyz = grass.read_command('v.out.ascii', input=vector_name, type='point', format='point', columns=rand_col, overwrite=True)\n",
    "    elev_list = [float(attr.split('|')[3]) if attr.split('|')[3] != '' else None for attr in xyz.split('\\n')[:-1]]\n",
    "    elev = np.asarray(elev_list, dtype=np.float64)\n",
    "    \n",
    "    return elev"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# aux func: fits a 4PL curve to mean of correlation values\n",
    "# plots and funcs from http://people.duke.edu/~ccc14/pcfb/analysis.html\n",
    "def logistic4(x, A, B, C, D):\n",
    "    ''' 4PL logistic equation ''' \n",
    "    return ((A-D)/(1.0+((x/C)**B))) + D\n",
    "\n",
    "def residuals(p, y, x):\n",
    "    ''' Deviations of data from fitted 4PL curve ''' \n",
    "    A,B,C,D = p\n",
    "    err = y-logistic4(x, A, B, C, D)\n",
    "    return err\n",
    "\n",
    "def peval(x, p):\n",
    "    ''' Evaluated value at x with current parameters ''' \n",
    "    A,B,C,D = p\n",
    "    return logistic4(x, A, B, C, D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# matplotlib figures appear inline in the notebook rather than in a new window.\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create GRASS GIS runtime environment\n",
    "# with this, you can run GRASS without startig a shell/gui session\n",
    "gisbase = subprocess.check_output([\"grass76\", \"--config\", \"path\"]).strip()\n",
    "os.environ['GISBASE'] = gisbase\n",
    "sys.path.append(os.path.join(gisbase, \"etc\", \"python\"))\n",
    "\n",
    "# GRASS GIS imports\n",
    "import grass.script as grass\n",
    "import grass.script.setup as gsetup\n",
    "import grass.script.array as garray\n",
    "import grass.script.vector as gvect"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set GRASS GIS session data\n",
    "# I use two systems, so this makes things a bit easier\n",
    "if sys.platform == \"linux\" or sys.platform == \"linux2\":\n",
    "    rcfile = gsetup.init(gisbase, \"/mnt/sda/grassdata/\", \"utm\", \"garopaba_22J\")\n",
    "elif sys.platform == \"darwin\":\n",
    "    rcfile = gsetup.init(gisbase, \"/Volumes/MacintoshHD2/grassdata/\", \"utm\", \"garopaba_22J\")\n",
    "# elif platform == \"win32\":\n",
    "    # Windows...\n",
    "    \n",
    "# grass.message('Current GRASS GIS 7 environment:')\n",
    "# print grass.gisenv()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# overwrite for GRASS modules\n",
    "ow = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data dir \n",
    "# use this to set different paths for different systems\n",
    "if sys.platform == \"linux\" or sys.platform == \"linux2\":\n",
    "    dataDir = '/mnt/sda/Dropbox/USP/projetosPesquisa/LiDAR_terrestre_SfM/_areas_estudo/garopaba/monteCarlo/'\n",
    "elif sys.platform == \"darwin\":\n",
    "    dataDir = '/Volumes/MacintoshHD2/Dropbox/USP/projetosPesquisa/LiDAR_terrestre_SfM/_areas_estudo/garopaba/monteCarlo/'\n",
    "#     dataDir = '_path_to_your_files_'\n",
    "os.chdir(dataDir)\n",
    "# os.getcwd()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# names for the files\n",
    "method='bilinear'\n",
    "step = 0.4\n",
    "dem_tls_10cm = 'tls_rinxyz_mean_10cm_' + method + '_step_' + str(step)\n",
    "dem_sfm_10cm = 'sfm_rinxyz_mean_10cm_' + method + '_step_' + str(step)\n",
    "diff_sfm_tls_10cm = 'diff_10cm_sfm_tls'\n",
    "mask_tls_sfm = 'mask_tls_sfm'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_random = 50 # n runs \n",
    "npoints_list = [50,100,250,500,1000,2500,5000,10000]\n",
    "dem_list = [dem_tls_10cm,dem_sfm_10cm]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run monte carlo sampling and save data as csv\n",
    "for dem,points in itertools.product(dem_list, npoints_list):\n",
    "    df = pd.DataFrame()\n",
    "    file_out = dem + '_rand_MC_' + str(points)\n",
    "    for run in range(n_random):\n",
    "        col_name = 'rand_' + str(points) + '_' + str(run).zfill(2)\n",
    "        elev = sample_dems_mc(dem, points, mc_run=run, ow_vector=True, ow_what=True, vmask=mask_tls_sfm)\n",
    "        df[col_name] = np.sort(elev)\n",
    "        df.to_csv(path_or_buf=file_out+'.csv', na_rep='NaN')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEICAYAAABfz4NwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xu8VXWd//HX+xwOiIYCwiiCiU5EoUOgWDZp3ipQywujpllmv2mcbjPORVSyUSQd8/KrXzY+apzS0kpFIqRGH9igjjaTJoaihCheCg4YeDmYcJLDOZ/fH+u7D5vNPpx14FzXeT8fj/U4a3+/a639/e61z/rs9f2u9V2KCMzMzGp6ugBmZtY7OCCYmRnggGBmZokDgpmZAQ4IZmaWOCCYmRnggNDtJI2VFJIGdOI2z5F0X2dtz3qepGMkrd5B/nck/Uu1ZSUtk3RMNxTTCsYBoRtIeknSh7pq+xHxo4j4SFdtv6gkHSvpAUkbJL1UJX9syt8k6Zmu3IcdFRGfi4ivtpF3cEQ82M1F2imSvi/pym54n+GSfippo6TfSfpEjnUGSlq+o8BcNA4IfVxnnmn0QxuBm4EZbeTfDiwB9gYuBeZKGtlNZetT+sD38EZgM7APcA7wbUkHt7PODGB9VxesV4kIT104AbcBLUAj8CZwERDAgJR/HvAC8EfgReCcdrZ3HvA/wDeAV4ErU9ovy5YJ4HPAc0AD2T+DUt63gZ+ULXsNsAgQ8DTwsbK8OuAVYDLwb6n8pWkLMKudsr5E9k+1lOzg+z2yf8h7U33/CxhWtvwRwP+mMj8JHFOW9xlgeVrvBeBvy/KOAVYD/wysA9YCn+nAPvoQ8FJF2juBt4AhZWkPA59L87XAl4HnU5keB/Yv+/y/kD7/PwJfBf481e0NYA4wsJ0yler05bQPXir/bgDfB64sX7bic/9Qmp+V3u/WVJZlwJSyZfcDfkJ24HsR+PuyvPcCv0r7Y236Dgwsyw/gi6meL+6gLiL7vq5L9X8KOAQ4H2giO1C/CfwsR5lmAXOBO1N9fgO8p53Pco/0Hu+s+L/82g7WOTB9304o/2yLPvV4AfrDVPEPOjb9Iw1IX9Q3gPEpbxRwcDvbOo/sYPx3aRuDqR4Qfg4MBd6e/rGmpbzdgWfTOkelg82YlHcRcGfZdk4BnqpShklpm5Nz1PsRsiAwOh0QfkMWYHYD7gcuT8uOJgtwJ5KduX44vR6Z8k8iO6gKOBrYBBya8o5Jn8lssiB2YsoftqPylZWzWkA4DVhekfZvwLfS/AyyA9v4VKb3AHuXff53A3sCB5MFlkXAQcBewG+BT7dTplKdvg4MSnXeWPZd+T75A8Kf0mdSC1wNPJLyasgC2WXAwFS+F4CpKf8wsiA9gOx7uxz4h4rv2S+A4cDgHdRlanqfoemzejcwqrIeOcs0iyyInJ729YVkQaNuB+8/GdhUkXYhKQC1sc7P03dgm8+26JObjHpeC3CIpMERsTYiluVYZ01EfCsitkREYxvLfC0iGiLi98ADZAdxImIT8CmyA80Pgb+LiFIb6Q+BEyXtmV5/iuyXVKvUZDI/rbckR1m/FRF/iIh6sl/Yj0bEkoj4E/BTsn9WgE8C90TEPRHREhG/ABaTHciIiP+MiOcj89/AfWQBraQJmB0RTRFxD9kvzvE5yteWtwEbKtI2AEPS/GeBr0TEilSmJyPi1bJlr42IN9L+fBq4LyJeiIgNZGdIk8nnXyLirVTn/wTO3Im6/DJ9rs1k+/M9Kf1wsoA7OyI2R8QLwH8AZwFExOMR8Uj6nr0E/DtZYCp3dUS8toPvIWT7ZgjwLrIz1eURsbaNZXdYpuTxiJgbEU1k3+PdyAJXW95G9sOrXPm+3Iak04DaiPjpDrZZSA4IPSgiNgIfJ2veWSvpPyW9K8eqq3Is83LZ/Cayf4rS+z5K9qtLZM0JpfQ1ZM1RfyVpKNnp8o9K+ZLqyE7XfxwRd+QoA8AfyuYbq7wulesA4AxJDaUJOJLsrAlJJ0h6RNJrKe9EYETZtl6NiC1t1XknvEn2C7/cnmTNFAD7kzUXtSVvvXfk9fQdKfkdWXNKR1V+F3ZLbf4HAPtVfOZfJjujQ9I7Jf1c0suS3gD+lW0/c8jxXYyI+8nOrm4E1km6qexHR6UdlqnyPSOihaxpbUefS3v7spWkPYBrgb/fca2KyQGhe7Q5pGxELIyID5Md+J4h+zW009vLQ9IXyZoh1pA1E5X7Admv9TOAX6Vf9iXfIvul9ZVdef82rAJui4ihZdMeEfE1SYPI2pSvB/aJiKHAPWQBrassAw6SVP4r8j0pvVTeP+/C9wcYlg5QJW8n22edZRVZ23/5Zz4kIk5M+d8m+06Oi4g9yQ7MlZ95ru9iRNwQEYcBE8j6Z0od+ZXrt1cmyIIxAJJqgDHs+HN5FhggaVxZWvm+LDeOrHnsYUkvA/OAUSkojm2nmn2eA0L3+ANZW+g2JO0j6ZT0T/8W2S+Zlq4siKR3knVEf5KsSegiSZPKFpkPHApcQNYRWVrvb8maC85Jv8o62w+Bj0maKqlW0m7p+voxZG3Jg8j6LbZIOgHY5ctsJdVI2o2sLVrpPQcCRMSzwBPA5Sn9NGAiWWAC+C7wVUnjlJkoae9dLVMVV6TLH48CPgrc1Ynb/jXwR0kXSxqcPvdDJB2e8oeQ/QB4M525fn5n3kTS4ZLel84wN5L1aZS+Q5X/G+2VCeAwSdPTWc4/kP3vPNLW+6ezrHnAbEl7SPoAWf/YbVUWf5os4ExK02dTGSeR78y8T3NA6B5XA19Jp7+nl6XXAP9E9uvmNbID7k790+WR/oF+CFyT2ryfI/vVd1v6FU5qC/4J2VUW88pWP5vsH3eNpDfT9OXOKltErCL7J/0y2YF/FdmvyJqI+CPZKfwc4HXgE8CCTnjbD5I139xD9uu7kaxvouQsYEp6z68Bp0dE6TLEr6fy3Ed20PweWQd/Z3o5vfcasqa7z0XEM5218dSn8FGyg92LZBcYfJes4xuyjtdPkDWt/AfZlT07Y8+0/utkzV6vAtelvO8BE1Lz0PwcZYKsw/7jaXufAqan/oQd+QLZ/llHdjnx50v9dZKOkvQmQOovebk0kf1ftqTXzTtZ/z6jdCmiWStJl5FdovfJni6LWTlJs4B3+LvZNXr7zSTWzSQNB/6a7JeXmfUjbjLqhZSNU/Nmlek7Xfy+f0PWVHNvRDyUY/m3t1HONyW9vSvLmpeycX2qle+cXlC2L7dRtnt7umwdVWp2qTZ10/v3+u9iX+AmIzMzA3yGYGZmSZ/qQxgxYkSMHTu2p4thZtanPP74469ERLsDM/apgDB27FgWL17c08UwM+tTJP0uz3JuMjIzM8ABwczMEgcEMzMDHBDMzCxxQDAzMyDnVUaSbiYbcGpdRBxSJV/AN9n6pKrzIuI3Ke/TbB0u+cqI+EFKP4zsaUmDyQYXuyB8l5z1IRuXrOONhS/R3PAWtUMHsefUsewx+c96ulhWILNmzcqV1lly3aks6YNkQzPf2kZAOJHskY4nAu8DvhkR70vj4iwmGzEyyB6Nd1hEvC7p12QjWD5KFhBuiIgd3rI/ZcqU8GWnXeeWf/48r63eOsLv8DH785n/++0eLFHvtXHJOl6/c8V26cM+Pt5BwTrFjg78HQ0Kkh6PiCntLZerySiNa/PaDhY5hSxYREQ8AgyVNIrsWaq/SI/Ye53s+avTUt6e6fF8QTbu/ql5ymJdozIYALy2ehW3/HOXjcbdp1ULBjtKN+sLOuvGtNFs+/CI1SltR+mrq6RvR9L5wPkAb3973xij6paLHmbTG1uHZ999zzo+c+1RO1ij51UGg/bSzax4en2nckTcFBFTImLKyJHt3nndadZecQXLDz6E5e96N8sPPoS1V1yRa73KYACw6Y0mbrno4a4opplZp+msgFBP2XNOyZ5xWt9O+pgq6b3C2iuuoOH2O6A5PSCpuZmG2+/IFRQqg0F76WZmvUVnBYQFwLnp2bJHABsiYi2wEPiIpGGShpE9B3dhyntD0hHpCqVzyR6L1ys03FH9SYFtpZuZFUHey05vB44BRkhaDVxO9mByIuI7ZFcJnQisJLvs9DMp7zVJXwUeS5uaHRGlzukvsPWy03vT1G2ufORK7nr2LlqihRrVcMY7z+ArR6SrY9u68spXxZpZN5k+fTrz5s2rmt5VcgWEiDi7nfwAvthG3s3AzVXSFwPbXcLaHa585EruXLH1135LtLS+bg0KZmY9aOLEiQAsWrSIDRs2sNdee3H88ce3pneFPjX8dWcpDwaV6Q4IlsfuR+zLpkderppu1lkmTpzYpQGgUr8MCGa7avip4wDY9OjL2S2Xgt3ft29rullf5IBgtpOGnzrOAcAKpdffh2BmZt3DAcEAGPS2IR1KN7PicUAwAI4/73xUW7tNmmprOf6883uoRGbW3dyHYAC8+6hjAXj4jlv546uvMGTvERx11rmt6WZWfA4I1urdRx3rAGDWj7nJyMzMgP50hrB0DiyaDRtWw9gxIPV0iczMepX+cYawdA7c/UXYsIrsLqKuM2b80A6lm5n1Fv0jINx7MTRv7pa3OuUfD93u4D9m/FBO+cdDu+X9zcx2Vv9oMmrc0dM/q5Cqj2yas5nJB38z64v6xxlCBw096+MdSjczK4L+cYagWojm3IuPuvxyABrm3JU9Na22lqFnntGabmZWRP0jIHQgGJSMuvxyBwAz61fcZGRmZoADgpmZJQ4IZmYG9NOA0NbFo2ozx8ys+PplQGjrXuXo4ruYzcx6s34ZEMzMbHuFDwjzl9TT0tOFMDPrAwodEOYvqWfmvKeQW4LMzNpV6IBw3cIVNDZ1/KY0M7P+qNABYU1DY08Xwcyszyh0QNhv6OCq6QOrjWQKDKwZ2JXFMTPr1QodEGZMHc/gutrt0me/8hqqCApCzP7A7O4qmplZr1Powe1OnTwagE13D2IP3tomrzaCLWXPN6jV9oHDzKw/KfQZAmRBYY/abc8GvjlsKFtqtq36ltjCN3/zze4smplZr5IrIEiaJmmFpJWSLqmSf4CkRZKWSnpQ0piyvGskPZ2mj5elf1/Si5KeSNOkzqlSFRWPz3x5QPWzgZc3vtxlRTAz6+3aDQiSaoEbgROACcDZkiZULHY9cGtETARmA1endU8CDgUmAe8DLpS0Z9l6MyJiUpqe2OXa5LTvluqXou67x77dVQQzs14nzxnCe4GVEfFCRGwG7gBOqVhmAnB/mn+gLH8C8FBEbImIjcBSYNquF3vXfHDTpqrPTP7gmA/2QGnMzHqHPAFhNLCq7PXqlFbuSWB6mj8NGCJp75Q+TdLukkYAxwL7l613VWpm+oakQTtVg3bMX1K/3ZB1D+2+O2j7kU0fWv1QVxTBzKxP6KxO5QuBoyUtAY4G6oHmiLgPuAf4X+B24FdAqb1mJvAu4HBgOHBxtQ1LOl/SYkmL169f36FCzV9Sz4y5T243vKn7EMzMtpcnINSz7a/6MSmtVUSsiYjpETEZuDSlNaS/V6U+gg+TPYrg2ZS+NjJvAbeQNU1tJyJuiogpETFl5MiRHarcFT9bRlPz9k1D7kMwM9tenoDwGDBO0oGSBgJnAQvKF5A0QlJpWzOBm1N6bWo6QtJEYCJwX3o9Kv0VcCrw9K5XZ1uvb2qqmn7B6w3s1rLtGKi71e7GBYde0NlFMDPrM9q9MS0itkj6ErAQqAVujohlkmYDiyNiAXAMcLWkAB4CvphWrwMezo75vAF8MiK2pLwfSRpJdtbwBPC5zqvWjp20cROQ3Y/wcl0d++6xLxccegEnHXRSdxXBzKzXyXWnckTcQ9YXUJ52Wdn8XGBulfX+RHalUbVtHtehku6EoYPraGhs4k9Ry2Bt20x00sZNnNS4BS7rWL+EmVlRFfpO5VknH0xdjWhU9UHuGPS27i2QmVkv1i/GMhp698bqCzS+3o2lMTPr3Qp9hlCyjhHVM/YaUz3dzKwfKnRAmL+knhl3Pcm/bj6DTVHxrIO6wXD8ZdVXNDPrhwodEGYtWEZTS7Cg5Ujuav4gW6KGCNhCDbznEzDxzJ4uoplZr1HogNDQmN2HcHLNLzm79kEGqAUJBtACS26DpXN6uIRmZr1HoQNCyay6WxmoLdsmNm+Ge6uOlmFm1i/1i4AwjDerZzS+1r0FMTPrxfpFQDAzs/b1i4CwkTZG1q7bo3sLYmbWixU6IGz/xAMzM2tLoQNCaeDrPXir+gJNbdzBbGbWDxU6INT4FMHMLLdCB4SW7Z+NY2ZmbSh0QCh5nTZGNR08vHsLYmbWixU6IAwdXAfArKZz2RIVVa2phROu6YFSmZn1ToUOCCOHbB3Qrnm7a44KXXUzsw4r9FHxuXXZVUQXDZjDoIonptHSBItm90CpzMx6p0IHhJL99Er1jA2rurcgZma9WL8ICC1tVVO13VsQM7NerF8EhFpaqmdEc/V0M7N+qF8EBF92ambWvn4REMI3qJmZtatfBIRhaut5CK93b0HMzHqxfhEQ1sSI6hl7jenegpiZ9WL9IiBcu+VM3oqKK4pq6uD4y3qmQGZmvVChA0Kttt6drMo7leWhUM3MyhU6IDSn3uSLBsxhoLZUZG72ncpmZmUKHRBKg9u1fafy6m4sjZlZ71bogFBqFXKnsplZ+wodEBo2NQGwqGVS9XsRxn2kewtkZtaL5QoIkqZJWiFppaRLquQfIGmRpKWSHpQ0pizvGklPp+njZekHSno0bfNOSQMrt7ur9kpNRh+rfaR6H/Kyn3b2W5qZ9VntBgRJtcCNwAnABOBsSRMqFrseuDUiJgKzgavTuicBhwKTgPcBF0raM61zDfCNiHgH8Drw17tencqyZ3+H0daNaa919luamfVZec4Q3gusjIgXImIzcAdwSsUyE4D70/wDZfkTgIciYktEbASWAtMkCTgOmJuW+wFw6s5Xo7pSk5GZmbUvT0AYDZQ/OGB1Siv3JDA9zZ8GDJG0d0qfJml3SSOAY4H9gb2BhojYsoNtAiDpfEmLJS1ev359njq12q2u0F0kZmadqrOOmBcCR0taAhwN1APNEXEfcA/wv8DtwK+ADo05HRE3RcSUiJgycuTIDhXqrS1tDHttZmbbyRMQ6sl+1ZeMSWmtImJNREyPiMnApSmtIf29KiImRcSHAQHPAq8CQyUNaGubnaHFo5yameWWJyA8BoxLVwUNBM4CFpQvIGmEpNK2ZgI3p/Ta1HSEpInAROC+iAiyvobT0zqfBu7e1cq0xc9DMDNrX7sBIbXzfwlYCCwH5kTEMkmzJZ2cFjsGWCHpWWAf4KqUXgc8LOm3wE3AJ8v6DS4G/knSSrI+he91Up2287PmI6rfh3DwaV31lmZmfc6A9heBiLiHrC+gPO2ysvm5bL1iqHyZP5FdaVRtmy+QXcHUZXavq2FTUwvH1zxR/T6E5+7ryrc3M+tTCn0ZTmNT1qnssYzMzNpX6IBQaiXyWEZmZu0rdEAoeSH2qd6HMPygbi+LmVlv1S8Cwl/WLK/eh/DSL7u9LGZmvVW/CAi1tHGDWnToHjkzs0IrdEDwQzLNzPIrdEAYOKDQ1TMz61SFPmKWxjKqb/Mqo/2rp5uZ9UOFDggl1245k01R8fydusFw/GXVVzAz64cKHRCGpiemLWg5kkuaPsvqlhG0hFjDCPjYDTDxzB4uoZlZ71HogDDr5IOpq8m6lhe0HMmRm29gfNOP+fUpDzkYmJlVyDWWUV916uTsmTvXLVzBmoZG9hs6mBlTx7emm5nZVoUOCJAFhVMnj4alc2DRxXD3anhwTNZ/4LMEM7NWhQ8IQBYMfvb30NSYvd6wKnsNDgpmZkmh+xBaLZq9NRiUNDVm6WZmBvSDgDB/ST0tbQ1z7eGvzcxaFTogzF9Sz8x5T7GmZe/qC3j4azOzVoUOCNctXEFjU7NvTDMzy6HQAWFNQ9ZvUHlj2uoW35hmZlap0FcZ7Td0MPVlQWHB5iMBGD10MP8z8bieLJqZWa9T6DOEGVPHU1e77SDYdbVixtTxPVQiM7Peq9ABAdj6YOW2XpuZGVDwgHDdwhU0tWwbAZpagusWruihEpmZ9V6FDgilTuW86WZm/Vm/6VQ+ueaXXDRgDvvpFdZpJCzd6KuMzMzKFPoMYcbU8dTViJNrfsnX6r7LmJpXqBHsy/psLKOlc3q6iGZmvUahAwIAgosGzGF3bd423WMZmZlto9AB4bqFK2hqDvbTK9UX8FhGZmatCh0QSp3Ha2JE9QU8lpGZWatCB4T9hg4G8FhGZmY55AoIkqZJWiFppaRLquQfIGmRpKWSHpQ0pizvWknLJC2XdIMkpfQH0zafSNOfdV61MjOmjmdwXe12YxltGjzKYxmZmVVoNyBIqgVuBE4AJgBnS5pQsdj1wK0RMRGYDVyd1v1L4APAROAQ4HDg6LL1zomISWlat6uVqXTq5NH81WGjqZVY0HIkRzd9i8smP8zuFz/jYGBmViHPGcJ7gZUR8UJEbAbuAE6pWGYCcH+af6AsP4DdgIHAIKAO+MOuFjqv+UvqufPXq2iO7G7l5gju/PUq5i+p764imJn1GXkCwmhgVdnr1Smt3JPA9DR/GjBE0t4R8SuyALE2TQsjYnnZerek5qJ/KTUldaZZC5ZVHbpi1oJlnf1WZmZ9Xmd1Kl8IHC1pCVmTUD3QLOkdwLuBMWRB5DhJR6V1zomIvwCOStOnqm1Y0vmSFktavH79+g4VqqGxqUPpZmb9WZ6AUA/sX/Z6TEprFRFrImJ6REwGLk1pDWRnC49ExJsR8SZwL/D+lF+f/v4R+DFZ09R2IuKmiJgSEVNGjhzZocqZmVl+eQLCY8A4SQdKGgicBSwoX0DSCEmlbc0Ebk7zvyc7cxggqY7s7GF5ej0irVsHfBR4eterY2ZmO6vdgBARW4AvAQuB5cCciFgmabakk9NixwArJD0L7ANcldLnAs8DT5H1MzwZET8j62BeKGkp8ATZGcd/dFqtkqGD6zqUbmbWnymi7zwxZsqUKbF48eLcy89fUs+Mu57cpmO5rkZcd8Z7OHVyZb+4mVkxSXo8Iqa0t1yhh78uHfSvW7iCNQ2N7Dd0MDOmjncwMDOrotABAbKg4ABgZta+Qo9lZGZm+TkgmJkZ4IBgZmaJA4KZmQEOCGZmljggmJkZ4IBgZmaJA4KZmQEOCGZmlhT+TuX5S+o9dIWZWQ6FDgjzl9Qzc95TNDY1A1Df0MjMeU8BOCiYmVUodJPRdQtXtAaDksamZq5buKKHSmRm1nsVOiCsaWjsULqZWX9W6ICw39DBHUo3M+vPCh0Qjn1X9Wcwt5VuZtafFTogPPDM+g6lm5n1Z4UOCO5DMDPLr9ABwX0IZmb5FTogzJg6nroabZNWVyNmTB3fQyUyM+u9Ch0QAFA7r83MDCh4QLhu4QqammObtKbm8I1pZmZVFDoguFPZzCy/QgcEdyqbmeVX6IAwY+p4BtfVbpM2uK7WncpmZlUUerTT0oimHv7azKx9hQ4IkAUFBwAzs/YVusnIzMzyc0AwMzPAAcHMzJJcAUHSNEkrJK2UdEmV/AMkLZK0VNKDksaU5V0raZmk5ZJukKSUfpikp9I2W9PNzKxntBsQJNUCNwInABOAsyVNqFjseuDWiJgIzAauTuv+JfABYCJwCHA4cHRa59vA3wDj0jRtVytjZmY7L88ZwnuBlRHxQkRsBu4ATqlYZgJwf5p/oCw/gN2AgcAgoA74g6RRwJ4R8UhEBHArcOou1cTMzHZJnoAwGlhV9np1Siv3JDA9zZ8GDJG0d0T8iixArE3TwohYntZf3c42AZB0vqTFkhavX+8H25iZdZXO6lS+EDha0hKyJqF6oFnSO4B3A2PIDvjHSTqqIxuOiJsiYkpETBk50o++NDPrKnluTKsH9i97PSaltYqINaQzBElvA/4qIhok/Q3wSES8mfLuBd4P3Ja20+Y2O8v8JfW+U9nMLIc8ZwiPAeMkHShpIHAWsKB8AUkjJJW2NRO4Oc3/nuzMYYCkOrKzh+URsRZ4Q9IR6eqic4G7O6E+25i/pJ6Z856ivqGRAOobGpk57ynmL+mS2GNm1qe1GxAiYgvwJWAhsByYExHLJM2WdHJa7BhghaRngX2Aq1L6XOB54CmyfoYnI+JnKe8LwHeBlWmZezulRmWuW7iCxqbmbdIam5r9PAQzsypyjWUUEfcA91SkXVY2P5fs4F+5XjPwt21sczHZpahdxs9DMDPLr9B3Kvt5CGZm+RU6IPh5CGZm+RV6+Gs/D8HMLL9CBwTw8xDMzPIqdJORmZnlV/gzBN+YZmaWT6EDQunGtNK9CKUb0wAHBTOzCoVuMvKNaWZm+RU6IPjGNDOz/AodEHxjmplZfoUOCL4xzcwsv0J3KvvGNDOz/AodEMA3ppmZ5VXoJiMzM8vPAcHMzAAHBDMzSxwQzMwMcEAwM7Ok8FcZeXA7M7N8Ch0QPLidmVl+hW4y8uB2Zmb5FTogeHA7M7P8Ch0QPLidmVl+hQ4IHtzOzCy/Qncqe3A7M7P8Ch0QwIPbmZnlVegmIzMzy88BwczMAAcEMzNLHBDMzAzIGRAkTZO0QtJKSZdUyT9A0iJJSyU9KGlMSj9W0hNl058knZryvi/pxbK8SZ1bNTMz64h2rzKSVAvcCHwYWA08JmlBRPy2bLHrgVsj4geSjgOuBj4VEQ8Ak9J2hgMrgfvK1psREXM7pyrVeXA7M7N88pwhvBdYGREvRMRm4A7glIplJgD3p/kHquQDnA7cGxGbdrawHVUa3K6+oZFg6+B285fUd1cRzMz6jDwBYTSwquz16pRW7klgepo/DRgiae+KZc4Cbq9Iuyo1M31D0qCcZc7Ng9uZmeXXWZ3KFwJHS1oCHA3UA61HYkmjgL8AFpatMxN4F3A4MBy4uNqGJZ0vabGkxevXr+9QoTy4nZlZfnkCQj2wf9nrMSmtVUSsiYjpETEZuDSlNZQtcibw04hoKltnbWTeAm4ha5raTkTcFBFTImLKyJEjc1WqxIPbmZnllycgPAaMk3SgpIFkTT8LyheQNEJSaVszgZsrtnE2Fc1F6awBSQJOBZ7uePF3zIPbmZnl1+5VRhGxRdKXyJp7aoGbI2KZpNnA4ohYABwDXC0gW8GRAAAE/ElEQVQpgIeAL5bWlzSW7Azjvys2/SNJIwEBTwCf2+XaVPDgdmZm+SkieroMuU2ZMiUWL17c08UwM+tTJD0eEVPaW853KpuZGeCAYGZmiQOCmZkBDghmZpY4IJiZGdDHrjKStB743U6sOgJ4pZOL09u5zv2D69w/7GqdD4iIdu/s7VMBYWdJWpznkqsicZ37B9e5f+iuOrvJyMzMAAcEMzNL+ktAuKmnC9ADXOf+wXXuH7qlzv2iD8HMzNrXX84QzMysHQ4IZmYGFDwgSJomaYWklZIu6eny7ApJ+0t6QNJvJS2TdEFKHy7pF5KeS3+HpXRJuiHVfamkQ8u29em0/HOSPt1TdcpLUq2kJZJ+nl4fKOnRVLc703M6kDQovV6Z8seWbWNmSl8haWrP1CQfSUMlzZX0jKTlkt5f9P0s6R/T9/ppSbdL2q2I+1nSzZLWSXq6LK3T9q2kwyQ9lda5IT1vJr+IKORE9uyG54GDgIFkz32e0NPl2oX6jAIOTfNDgGeBCcC1wCUp/RLgmjR/InAv2fMmjgAeTenDgRfS32FpflhP16+duv8T8GPg5+n1HOCsNP8d4PNp/gvAd9L8WcCdaX5C2v+DgAPT96K2p+u1g/r+APhsmh8IDC3yfiZ7RvuLwOCy/XteEfcz8EHgUODpsrRO27fAr9OySuue0KHy9fQH1IUf/PuBhWWvZwIze7pcnVi/u4EPAyuAUSltFLAizf87cHbZ8itS/tnAv5elb7Ncb5vIHtm6CDgO+Hn6or8CDKjcz2QPcXp/mh+QllPlvi9frrdNwF7p4KiK9MLu5xQQVqUD3IC0n6cWdT8DYysCQqfs25T3TFn6NsvlmYrcZFT6kpWsTml9XjpFngw8CuwTEWtT1svAPmm+rfr3tc/l/wEXAS3p9d5AQ0RsSa/Ly99at5S/IS3fl+p8ILAeuCU1k31X0h4UeD9HRD1wPfB7YC3ZfnucYu/ncp21b0en+cr03IocEApJ0tuAnwD/EBFvlOdF9rOgMNcRS/oosC4iHu/psnSjAWRNCt+OiMnARrJmhFYF3M/DgFPIguF+wB7AtB4tVA/p6X1b5IBQT/Ys55IxKa3PklRHFgx+FBHzUvIfJI1K+aOAdSm9rfr3pc/lA8DJkl4C7iBrNvomMFRS6Xng5eVvrVvK3wt4lb5V59XA6oh4NL2eSxYgiryfPwS8GBHrI6IJmEe274u8n8t11r6tT/OV6bkVOSA8BoxLVyoMJOt8WtDDZdpp6WqB7wHLI+LrZVkLgNJVBp8m61sopZ+brlQ4AtiQTksXAh+RNCz9MvtISut1ImJmRIyJiLFk++/+iDgHeAA4PS1WWefSZ3F6Wj5S+lnp6pQDgXFknW+9TkS8DKySND4lHQ/8lgLvZ7KmoiMk7Z6+56U6F3Y/V+iUfZvy3pB0RPoczy3bVj493cHSxZ03J5JdjfM8cGlPl2cX63Ik2ankUuCJNJ1I1na6CHgO+C9geFpewI2p7k8BU8q29X+AlWn6TE/XLWf9j2HrVUYHkf2jrwTuAgal9N3S65Up/6Cy9S9Nn8UKOnjlRQ/UdRKwOO3r+WRXkhR6PwNXAM8ATwO3kV0pVLj9DNxO1k/SRHY2+NeduW+BKekzfB74NyouTmhv8tAVZmYGFLvJyMzMOsABwczMAAcEMzNLHBDMzAxwQDAzs8QBwczMAAcEMzNL/j/FPozttr25qAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x104278ad0>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# reads data from csv files and calculates correlation\n",
    "dem = dem_tls_10cm\n",
    "avg_corr = []\n",
    "df_corr = pd.DataFrame()\n",
    "for points in npoints_list:\n",
    "    csv_file = dem + '_rand_MC_' + str(points) + '.csv'\n",
    "    #\n",
    "    df = pd.read_csv(csv_file, index_col=0)\n",
    "    # correlation of first column[0] with  all the others [1:]. \n",
    "    # No need to define column by name\n",
    "    corr = df.corr().iloc[0,1:]\n",
    "    avg_corr.append(corr.mean())\n",
    "    #\n",
    "    # plot correlation values for this set of random points\n",
    "    x_ax = np.empty(n_random -1)\n",
    "    x_ax.fill(points)\n",
    "    plt.plot(x_ax, corr, 'o')\n",
    "    plt.title(dem)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,u'tls_rinxyz_mean_10cm_bilinear_step_0.4')"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEICAYAAABfz4NwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3XucXVV99/HPdzIzCZEkkBCSdIabnYCNglwCYq1yebwAtSBIK0gUsC1Say9PH1SQllqUoj6+6vNgeWkpInKpQKMo0vAKlpv6FJBQ7obEIQGSYSJDAgkhwNx+zx97nXHnMHPOSTKTc2b29/167dfZZ+3LWevsc9Zvr7X22UcRgZmZWVO9M2BmZo3BAcHMzAAHBDMzSxwQzMwMcEAwM7PEAcHMzAAHhJ1O0r6SQlLzKO7zDEm3j9b+rP4kHS1pbYXl35L0d8OtK+kJSUfvhGzaBOOAsBNIelrSe8dq/xFxfUS8f6z2P1FJOkbSXZI2Snp6mOX7puVbJD05lsdwW0XEuRHxxRGWvTUi7t7JWdoukq6W9KWd8DozJd0s6RVJz0j6aA3btEpaXikwTzQOCOPcaLY0CugV4CrgMyMs/x7wEDALuBBYLGn2TsrbuDIOPoeXA73AHOAM4JuS3lplm88APWOdsYYSEZ7GcAKuBQaBV4HNwGeBAJrT8rOAVcDLwGrgjCr7Owv4f8DXgfXAl1Laz3PrBHAu8CvgJbIvg9KybwLfz637FeAOQMDjwB/klrUALwCHAP+c8l+a+oEvVMnr02RfqkfJKt9vk30hb0vl/U9g99z6RwL/lfL8CHB0btnZwPK03Srgk7llRwNrgf8FPA90A2dvwzF6L/B0Wdr+wOvAtFzaz4Bz0/wk4PPAUylPDwJ75d7/T6X3/2Xgi8Bvp7JtAm4CWqvkqVSmz6dj8HT+swFcDXwpv27Z+/7eNP+F9HrXpLw8ASzMrftbwPfJKr7VwF/mlh0B3JuOR3f6DLTmlgfw56mcqyuURWSf1+dT+R8D3gacA/SRVdSbgR/XkKcvAIuBG1N5/ht4e5X38k3pNfYv+15+ucI2+6XP2/H593aiT3XPQBGmsi/ovumL1Jw+qJuAA9KyecBbq+zrLLLK+C/SPnZh+IBwK7AbsHf6Yh2Xlk0FVqZt3p0qm/a07LPAjbn9nAQ8NkweDk77PKSGct9HFgTaUoXw32QBZgpwJ/D3ad02sgB3AlnL9X3p+ey0/PfJKlUBRwFbgEPTsqPTe3IxWRA7IS3fvVL+cvkcLiCcDCwvS/tn4Btp/jNkFdsBKU9vB2bl3v8fAdOBt5IFljuANwMzgF8CZ1bJU6lM/wRMTmV+JfdZuZraA8Jr6T2ZBFwK3JeWNZEFsouA1pS/VcAH0vLDyIJ0M9nndjnw12Wfs58AM4FdKpTlA+l1dkvv1e8A88rLUWOevkAWRE5Nx/o8sqDRUuH1DwG2lKWdRwpAI2xza/oMbPXeTvTJXUb1Nwi8TdIuEdEdEU/UsM1zEfGNiOiPiFdHWOfLEfFSRDwL3EVWiRMRW4CPkVU01wF/ERGlPtLrgBMkTU/PP0Z2JjUkdZn8MG33UA15/UZE/DoiusjOsO+PiIci4jXgZrIvK8AiYElELImIwYj4CbCMrCIjIv4jIp6KzD3A7WQBraQPuDgi+iJiCdkZ5wE15G8kuwIby9I2AtPS/J8AfxsRK1KeHomI9bl1vxoRm9LxfBy4PSJWRcRGshbSIdTm7yLi9VTm/wD+aDvK8vP0vg6QHc+3p/TDyQLuxRHRGxGrgH8FTgOIiAcj4r70OXsa+BeywJR3aURsqPA5hOzYTAPeQtZSXR4R3SOsWzFPyYMRsTgi+sg+x1PIAtdIdiU78crLH8utSDoZmBQRN1fY54TkgFBHEfEK8BGy7p1uSf8h6S01bLqmhnXW5ea3kH0pSq97P9lZl8i6E0rpz5F1R31Y0m5kzeXrS8sltZA11/8tIm6oIQ8Av87NvzrM81K+9gH+UNJLpQn4PbJWE5KOl3SfpA1p2QnAHrl9rY+I/pHKvB02k53h500n66YA2Iusu2gktZa7khfTZ6TkGbLulG1V/lmYkvr89wF+q+w9/zxZiw5J+0u6VdI6SZuAf2Tr9xxq+CxGxJ1kravLgeclXZE76ShXMU/lrxkRg2Rda5Xel2rHcoikNwFfBf6ycqkmJgeEnWPEW8pGxNKIeB9Zxfck2dnQdu+vFpL+nKwb4jmybqK875Kdrf8hcG86sy/5BtmZ1t/uyOuPYA1wbUTslpveFBFfljSZrE/5a8CciNgNWEIW0MbKE8CbJeXPIt+e0kv5/e0xfH2A3VMFVbI32TEbLWvI+v7z7/m0iDghLf8m2WdyfkRMJ6uYy9/zmj6LEXFZRBwGLCAbnykN5JdvXy1PkAVjACQ1Ae1Ufl9WAs2S5ufS8scybz5Z99jPJK0DfgDMS0Fx3yrFHPccEHaOX5P1hW5F0hxJJ6Uv/etkZzKDY5kRSfuTDUQvIusS+qykg3Or/BA4FPgrsoHI0nafJOsuOCOdlY2264A/kPQBSZMkTUnX17eT9SVPJhu36Jd0PLDDl9lKapI0hawvWuk1WwEiYiXwMPD3Kf1k4CCywARwJfBFSfOVOUjSrB3N0zD+IV3++G7gg8C/j+K+fwG8LOlzknZJ7/vbJB2elk8jOwHYnFquf7Y9LyLpcEnvSC3MV8jGNEqfofLvRrU8ARwm6ZTUyvlrsu/OfSO9fmpl/QC4WNKbJL2LbHzs2mFWf5ws4Bycpj9JeTyY2lrm45oDws5xKfC3qfl7ai69CfgbsrObDWQV7nZ96WqRvkDXAV9Jfd6/IjvruzadhZP6gr9PdpXFD3Kbn072xX1O0uY0fX608hYRa8i+pJ8nq/jXkJ1FNkXEy2RN+JuAF4GPAreMwsu+h6z7ZgnZ2ferZGMTJacBC9Nrfhk4NSJKlyH+U8rP7WSV5rfJBvhH07r02s+Rdd2dGxFPjtbO05jCB8kqu9VkFxhcSTbwDdnA60fJulb+lezKnu0xPW3/Ilm313rgf6dl3wYWpO6hH9aQJ8gG7D+S9vcx4JQ0nlDJp8iOz/NklxP/WWm8TtK7JW0GSOMl60oT2fdyMD0f2M7yjxulSxHNhki6iOwSvUX1zotZnqQvAB3+bI6NRv8xie1kkmYCf0x25mVmBeIuowak7D41m4eZvjXGr/unZF01t0XET2tYf+8R8rlZ0t5jmddaKbuvz3D5O6MB8vb5EfJ2W73ztq1K3S7DTTvp9Rv+szgeuMvIzMwAtxDMzCwZV2MIe+yxR+y77771zoaZ2bjy4IMPvhARVW/MOK4Cwr777suyZcvqnQ0zs3FF0jO1rOcuIzMzAxwQzMwscUAwMzPAAcHMzBIHBDMzA2oMCJKukvS8pMdHWC5Jl0nqlPSopENzy86U9Ks0nZlLP0zSY2mbyySN5a2MzczqanBwkJUrV3LPPfewcuVKBgfH9MbG26XWy06vJvuDi2tGWH482X3E5wPvILuP+jvSfXH+nuyOkQE8KOmWiHgxrfOnwP1kd5s8juyfpMwayuDgIJ2dnXR3dzNv3jw6OjpoanLj2mo3ODjIddddR1dXF729vbS2ttLW1saiRYsa6rNUU0CIiJ+q8p9DnARcE9l9MO6TtJukeWT/R/qTiNgAIOknwHGS7gamR8R9Kf0a4EM4IFiDGS9fZGtsnZ2dQ58hgN7eXrq6uujs7GT//fevc+5+Y7Q+0W1s/ecRa1NapfS1w6S/gaRzJC2TtKynp2e4VWwHjYembL1U+iKb1aq7u3voM1TS29vLunXrRtiiPhr+l8oRcQVwBcDChQt9J75RFBH09/dz/fXX89xzz9HX10dLSwtz587lxBNPBGBgYIDBwcGhx4ggIobmh0sbzfnSzRcrzVdbnr+B47bus6enZ9gv8tKlS/nFL35R8bXy+xqN52O573q+VhH09/cPm37vvffywAMPVN1+99135xOf+MRoZ+sNRisgdJH7n1Oy/zjtStPRZel3p/T2YdYvtIhgYGCA119/fWjq7e19w3xfXx/9/f1Dj8NN5ctGOusfGBjg1VdfHXre19fH2rVrufHGG9l1111pampi0qRJNDU1DU2SkLTVfPnzWpeV9ll6jXw6sNU+8s9L88Ol1TJf67pr1qzhrrvu2uoL3dzczBFHHME+++wz4nb516nl+basO5b7qmc+JrLx0vU4WgHhFuDTkm4gG1TeGBHdkpYC/yhp97Te+4ELImKDpE2SjiQbVP442R+4T0ivv/46mzZtYtOmTWzcuJFNmzaxZcsWXn31VbZs2TI0DQwM0NzczJQpU2htbWXy5MlMnjx5aL70OHXqVFpaWmhubt5qKk8rPc9XtuXuuece7r777q3SIoIDDzyQ97znPTvpHWpcbW1tW3Ublb7Ihx9+eEN9ka2xNTU1sWjRIjo7O1m3bh1z585tyIsTagoIkr5Hdqa/h6S1ZFcOtQBExLfIrhI6AegEtgBnp2UbJH0RKLWJLi4NMJP9x+nVZP9zehsTYEB506ZNdHd309PTw6pVq3jhhRfo6+tjypQpzJgxgxkzZjBt2jSmT5/OnnvuydSpU4emXXbZhUmTJu30PM+bN4/W1tatukVaW1uZO3fuTs9LIxovX2RrfE1NTey///4NNYhcblz9Qc7ChQujke52umHDBlauXMmaNWvo6upiypQpzJ07l2effZaXX36Z/v5+WlpaaG9vb7imYcl4acqa2faT9GBELKy2XsMPKjeaF198kYceeognn3yS5uZm5s+fz2GHHcaJJ57I5MmTWblyJcuXLx/qc+7r62vIy8tKfAZsZiUOCBXkf5DU0tLCqlWr2LhxI4cccghnnHEGM2bMeMM2lS4va8SAAOOjKWtmY88BYQSlrpS1a9fS19cHwJw5czj33HMr9vW7T97Mxiv3C4ygs7OTZ599digYQNZd9NRTT1XcrqOjg7a2NlpbWwGG+uQ7OjrGNL9mZjvKLYRhRAT33XcfAwMDW6XX0vXjPnkzG68cEIbx85//nC1bttDS0rJVC6HWrh/3yZvZeOTT1jKPPfYYTz75JGeddRbt7e3u+jGzwnALIWfjxo3ceeednH322UyZMsVdP2ZWKA4ISUTw4x//mGOOOYbp06cD7voxs2IpdEDI/86gqamJvr4+DjzwwHpny8ysLgobEMpv2SCJOXPmEBGFugujmVlJYTvEy//4JCJYv369//jEzAqrsAFhuFtM9PX1Ndw/GJmZ7SyFDQilW0zk+RYTZlZkhQ0IpVtMlO5L5N8ZmFnRFXZQuXSLiWuvvZZdd92VAw880L8zMLNCK2xAgOw/XTds2MDpp5/+hu4jM7OiKfTp8Pr165kxY4aDgZkZBQ8IzzzzDHvvvXe9s2Fm1hAKHRDWrFnjgGBmlhQ6IKxbt4558+bVOxtmZg2hsAGhv7+fLVu2sOuuu9Y7K2ZmDaGwAaGnp4c999zT9y0yM0sKGxCef/559txzz3pnw8ysYRQ2IGzYsIGZM2fWOxtmZg2j0AFh1qxZ9c6GmVnDKGxAWL9+vVsIZmY5hQwIEcGmTZuG/irTzMwKGhBee+01dtllF19hZGaWU8iA4NaBmdkbFTIgvPzyy0ybNq3e2TAzayiFDAibNm1yQDAzK1NTQJB0nKQVkjolnT/M8n0k3SHpUUl3S2rPLfuKpMfT9JFc+tWSVkt6OE0Hj06RqnMLwczsjaoGBEmTgMuB44EFwOmSFpSt9jXgmog4CLgYuDRt+/vAocDBwDuA8yTlO+8/ExEHp+nhHS5NjV5++WWPIZiZlamlhXAE0BkRqyKiF7gBOKlsnQXAnWn+rtzyBcBPI6I/Il4BHgWO2/Fs7xi3EMzM3qiWgNAGrMk9X5vS8h4BTknzJwPTJM1K6cdJmippD+AYYK/cdpekbqavS5o83ItLOkfSMknLenp6ashudZs3b/ZdTs3MyozWoPJ5wFGSHgKOArqAgYi4HVgC/BfwPeBeYCBtcwHwFuBwYCbwueF2HBFXRMTCiFg4e/bsUcnsli1bmDp16qjsy8xsoqglIHSx9Vl9e0obEhHPRcQpEXEIcGFKeyk9XpLGCN4HCFiZ0rsj8zrwHbKuqZ2ir6+PlpaWnfVyZmbjQi0B4QFgvqT9JLUCpwG35FeQtIek0r4uAK5K6ZNS1xGSDgIOAm5Pz+elRwEfAh7f8eJUNzAwwKRJk3bGS5mZjSvN1VaIiH5JnwaWApOAqyLiCUkXA8si4hbgaOBSSQH8FPjztHkL8LN0i4hNwKKI6E/Lrpc0m6zV8DBw7ugVa2Svvvoqu+yyy854KTOzcaVqQACIiCVkYwH5tIty84uBxcNs9xrZlUbD7fPYbcrpKHFAMDMbXuF+qeyAYGY2vMIFhC1btjggmJkNo1ABYXBwkNWrV/PCCy+wcuVKBgcH650lM7OGUdMYwkQwODjIddddx7PPPsvAwADr1q2jra2NRYsW0dRUqLhoZjaswtSEnZ2ddHV1MTCQ/S6ut7eXrq4uOjs765wzM7PGUJiA0N3dTW9v71Zpvb29rFu3rk45MjNrLIUJCPPmzaO1tXWrtNbWVubOnVunHJmZNZbCBISOjg7a2tqG/ke5tbWVtrY2Ojo66pwzM7PGUJhB5aamJhYtWsSVV15Je3s7HR0ddHR0eEDZzCwpTECALCi0tLRw5JFHMnPmzHpnx8ysoRTu9Pj1119/w1iCmZkVMCD09vY6IJiZDaNwAcH/hWBmNrzCBQRg6EojMzP7jUIFhMHBQQcDM7MRFCog9PX1efzAzGwEhQoIvsLIzGxkhQoIvsLIzGxkhQsIkydPrnc2zMwaUqECgi85NTMbWaECQn9/P83Nhbpbh5lZzQoXECZNmlTvbJiZNaTCBQS3EMzMhueAYGZmQMECwsDAgAOCmdkIChUQPIZgZjaywgUEtxDMzIbngGBmZoADgpmZJQ4IZmYGFCwg+CojM7ORFSog+CojM7OR1RQQJB0naYWkTknnD7N8H0l3SHpU0t2S2nPLviLp8TR9JJe+n6T70z5vlDTm96V2l5GZ2ciqBgRJk4DLgeOBBcDpkhaUrfY14JqIOAi4GLg0bfv7wKHAwcA7gPMkTU/bfAX4ekR0AC8Cf7zjxanMXUZmZiOrpYVwBNAZEasiohe4ATipbJ0FwJ1p/q7c8gXATyOiPyJeAR4FjlP2x8bHAovTet8FPrT9xaiNWwhmZiOrJSC0AWtyz9emtLxHgFPS/MnANEmzUvpxkqZK2gM4BtgLmAW8FBH9FfYJgKRzJC2TtKynp6eWMo3IAcHMbGSjNah8HnCUpIeAo4AuYCAibgeWAP8FfA+4FxjYlh1HxBURsTAiFs6ePXuHMumAYGY2sloCQhfZWX1Je0obEhHPRcQpEXEIcGFKeyk9XhIRB0fE+wABK4H1wG6Smkfa51hwQDAzG1ktAeEBYH66KqgVOA24Jb+CpD0klfZ1AXBVSp+Uuo6QdBBwEHB7RATZWMOpaZszgR/taGGq8WWnZmYjqxoQUj//p4GlwHLgpoh4QtLFkk5Mqx0NrJC0EpgDXJLSW4CfSfolcAWwKDdu8DngbyR1ko0pfHuUyjQiX2VkZjaymmrHiFhCNhaQT7soN7+Y31wxlF/nNbIrjYbb5yqyK5jG3ODgIJ2dnWzevJnVq1dzwAEH0NRUqN/kmZlVNeFPlwcHB7nuuuvo6uqit7eXm2++mfb2dhYtWuSgYGaWM+FrxM7OzqFgANDX10dXVxednZ11zpmZWWOZ8AGhu7t7KBiU9Pb2sm7dujrlyMysMU34gDBv3jxaW7e+TVJraytz586tU47MzBrThA8IHR0dtLW1DQWF1tZW2tra6OjoqHPOzMway4QfVG5qamLRokWsXLmSm2++mQ9/+MN0dHR4QNnMrMyEDwiQBYWOjg6mTp3K/vvvX+/smJk1pMKcJg8ODrpVYGZWQWFqSAcEM7PKClNDOiCYmVVWmBpyYGDAAcHMrILC1JCDg4O+06mZWQWFCghuIZiZjawwNaQDgplZZYWpIR0QzMwqK0wN6UFlM7PKClNDelDZzKyyQgUEtxDMzEZWmBrSAcHMrLLC1JAOCGZmlRWmhvSgsplZZYWpId1CMDOrrDA1pAOCmVllhakhHRDMzCorTA3p3yGYmVVWmIDgQWUzs8oKU0O6y8jMrLLC1JAOCGZmlRWmhnRAMDOrrDA1pAeVzcwqK1RAcAvBzGxkhakhfZWRmVllNdWQko6TtEJSp6Tzh1m+j6Q7JD0q6W5J7bllX5X0hKTlki6TpJR+d9rnw2nac/SK9UZuIZiZVVa1hpQ0CbgcOB5YAJwuaUHZal8DromIg4CLgUvTtr8LvAs4CHgbcDhwVG67MyLi4DQ9v6OFqcQBwcysslpqyCOAzohYFRG9wA3ASWXrLADuTPN35ZYHMAVoBSYDLcCvdzTT28ODymZmldUSENqANbnna1Na3iPAKWn+ZGCapFkRcS9ZgOhO09KIWJ7b7jupu+jvSl1J5SSdI2mZpGU9PT01ZHd4biGYmVU2WjXkecBRkh4i6xLqAgYkdQC/A7STBZFjJb07bXNGRBwIvDtNHxtuxxFxRUQsjIiFs2fP3u4MelDZzKyyWmrILmCv3PP2lDYkIp6LiFMi4hDgwpT2Ellr4b6I2BwRm4HbgHem5V3p8WXg38i6psaMWwhmZpXVUkM+AMyXtJ+kVuA04Jb8CpL2kFTa1wXAVWn+WbKWQ7OkFrLWw/L0fI+0bQvwQeDxHS/OyBwQzMwqq1pDRkQ/8GlgKbAcuCkinpB0saQT02pHAyskrQTmAJek9MXAU8BjZOMMj0TEj8kGmJdKehR4mKzF8a+jVqpheFDZzKyy5lpWioglwJKytIty84vJKv/y7QaATw6T/gpw2LZmdke4hWBmVllhakgHBDOzygpTQzogmJlVVpga0pedmplVVpga0i0EM7PKClNDRoQDgplZBYWpIQcHBxnh7hhmZkaBAoJbCGZmlRWmhnQLwcysssIEBLcQzMwqK0wN6RaCmVllhQkIEeGAYGZWQaECgruMzMxGVpga0l1GZmaVFSYguIVgZlZZYWpItxDMzCorTEBwC8HMrLLC1JC+ysjMrLLCBAR3GZmZVVaYgOAuIzOzygpTQ7qFYGZWWWECglsIZmaVFaaGdAvBzKyywgQEtxDMzCorTA0ZEfXOgplZQytMQADcZWRmVkGhAoKZmY3MAcHMzAAHBDMzSxwQzMwMcEAwM7PEAcHMzAAHBDMzS2oKCJKOk7RCUqek84dZvo+kOyQ9KuluSe25ZV+V9ISk5ZIuU/oxgKTDJD2W9jmUPhb8ozQzs+qqBgRJk4DLgeOBBcDpkhaUrfY14JqIOAi4GLg0bfu7wLuAg4C3AYcDR6Vtvgn8KTA/TcftaGFG4ttWmJlVV0steQTQGRGrIqIXuAE4qWydBcCdaf6u3PIApgCtwGSgBfi1pHnA9Ii4L7LT92uAD+1QSSrwje3MzKqrJSC0AWtyz9emtLxHgFPS/MnANEmzIuJesgDRnaalEbE8bb+2yj5HjVsIZmbVjVYteR5wlKSHyLqEuoABSR3A7wDtZBX+sZLevS07lnSOpGWSlvX09GxX5vx/ymZm1dUSELqAvXLP21PakIh4LiJOiYhDgAtT2ktkrYX7ImJzRGwGbgPembZvr7TP3L6viIiFEbFw9uzZNRZra+4yMjOrrpaA8AAwX9J+klqB04Bb8itI2kNSaV8XAFel+WfJWg7NklrIWg/LI6Ib2CTpyHR10ceBH41CeYblLiMzs+qq1pIR0Q98GlgKLAduiognJF0s6cS02tHACkkrgTnAJSl9MfAU8BjZOMMjEfHjtOxTwJVAZ1rntlEp0TDcQjAzq665lpUiYgmwpCztotz8YrLKv3y7AeCTI+xzGdmlqGPOLQQzs+oKUUu6hWBmVl0hAoJbCGZm1RWilvRlp2Zm1RUiILjLyMysukIEBHcZmZlVV4ha0i0EM7PqChEQ3EIwM6uuELWkWwhmZtUVIiC4hWBmVl0haklfdmpmVl0hAoK7jMzMqitEQHCXkZlZdYWoJd1CMDOrrhABwS0EM7PqClFLelDZzKy6QgQEdxmZmVVXiIDgLiMzs+oKUUu6hWBmVl0hAoLHEMzMqitEQBgcHHSXkZlZFYWoJd1CMDOrrjABwS0EM7PKClFLelDZzKy6QgQEdxmZmVVXiIDgQWUzs+oKUUu6hWBmVl0hAoJbCGZm1RWilnQLwcysuuZ6Z2Bn2G+//ejr66t3NszMGlohAsKsWbPqnQUzs4ZXiC4jMzOrzgHBzMwABwQzM0tqCgiSjpO0QlKnpPOHWb6PpDskPSrpbkntKf0YSQ/nptckfSgtu1rS6tyyg0e3aGZmti2qDipLmgRcDrwPWAs8IOmWiPhlbrWvAddExHclHQtcCnwsIu4CDk77mQl0ArfntvtMRCwenaKYmdmOqKWFcATQGRGrIqIXuAE4qWydBcCdaf6uYZYDnArcFhFbtjezZmY2dmoJCG3AmtzztSkt7xHglDR/MjBNUvm1nqcB3ytLuyR1M31d0uThXlzSOZKWSVrW09NTQ3bNzGx7jNag8nnAUZIeAo4CuoCB0kJJ84ADgaW5bS4A3gIcDswEPjfcjiPiiohYGBELZ8+ePUrZNTOzcrX8MK0L2Cv3vD2lDYmI50gtBEm7Ah+OiJdyq/wRcHNE9OW26U6zr0v6DllQqejBBx98QdIzNeR5OHsAL2zntuOJyzmxuJwTTz3Kuk8tK9USEB4A5kvajywQnAZ8NL+CpD2ADRExSHbmf1XZPk5P6flt5kVEt7KbDH0IeLxaRiJiu5sIkpZFxMLt3X68cDknFpdz4mnkslbtMoqIfuDTZN09y4GbIuIJSRdLOjGtdjSwQtJKYA5wSWl7SfuStTDuKdv19ZIeAx4ji5hf2qGSmJnZDqnpXkYRsQRYUpZ2UW5+MTDs5aMR8TRvHIQmIo7dloyamdnYKtIvla+odwZ2EpdzYnE5J56GLasiot55MDOzBlCkFoKZmVXggGBmZkBBAkK1m/M1Mkl7SbpL0i8lPSHpr1L6TEk/kfSr9Lh7Speky1JZH5V0aG5fZ6b6kmxqAAAECklEQVT1fyXpzHqVqRJJkyQ9JOnW9Hw/Sfen8twoqTWlT07PO9PyfXP7uCClr5D0gfqUpDJJu0laLOlJScslvXMiHlNJ/zN9bh+X9D1JUybCMZV0laTnJT2eSxu14yfpMEmPpW0uS5fnj72ImNATMAl4Cngz0Ep2m40F9c7XNuR/HnBomp8GrCS7d9RXgfNT+vnAV9L8CcBtgIAjgftT+kxgVXrcPc3vXu/yDVPevwH+Dbg1Pb8JOC3Nfwv4szT/KeBbaf404MY0vyAd48nAfunYT6p3uYYp53eBP0nzrcBuE+2Ykl1duBrYJXcsz5oIxxR4D3Ao8HgubdSOH/CLtK7StsfvlHLV+0OzEw7cO4GluecXABfUO187UJ4fkd15dgUwL6XNA1ak+X8BTs+tvyItPx34l1z6Vus1wkT2K/g7gGOBW9OX4QWgufxYkv0u5p1pvjmtp/Ljm1+vUSZgRqooVZY+oY4pv7kP2sx0jG4FPjBRjimwb1lAGJXjl5Y9mUvfar2xnIrQZVTLzfnGhdSEPgS4H5gTv7n9xzqyHwTCyOUdD+/D/wE+Cwym57OAlyL7cSRsneeh8qTlG9P646Gc+wE9wHdS99iVkt7EBDumEdFFdmv8Z4FusmP0IBPzmMLoHb+2NF+ePuaKEBAmBGX3iPo+8NcRsSm/LLLTiHF9/bCkDwLPR8SD9c7LTtBM1t3wzYg4BHiFrIthyAQ5pruT3Qp/P+C3gDcBx9U1UzvJeD1+RQgIVW/O1+gktZAFg+sj4gcp+dfK7iJbupvs8yl9pPI2+vvwLuBESU+T/efGscD/BXaTVPpFfT7PQ+VJy2cA62n8ckJ2xrc2Iu5PzxeTBYiJdkzfC6yOiJ7Ibmz5A7LjPBGPKYze8etK8+XpY64IAWHo5nzpaobTgFvqnKeapasLvg0sj4h/yi26BShdlXAm2dhCKf3j6cqGI4GNqRm7FHi/pN3Tmdv72fp25HUVERdERHtE7Et2jO6MiDPI/nDp1LRaeTlL5T81rR8p/bR0xcp+wHyyAbqGERHrgDWSDkhJ/wP4JRPsmJJ1FR0paWr6HJfKOeGOaTIqxy8t2yTpyPS+fTy3r7FV74GZnTT4cwLZ1TlPARfWOz/bmPffI2t6Pgo8nKYTyPpW7wB+BfwnMDOtL7K/PH2K7MaBC3P7+gTZ35h2AmfXu2wVynw0v7nK6M1kX/5O4N+BySl9SnremZa/Obf9han8K9hJV2dsRxkPBpal4/pDsqtMJtwxBf4BeJLsbsbXkl0pNO6PKdmffXUDfWQtvj8ezeMHLEzv2VPAP1N2AcJYTb51hZmZAcXoMjIzsxo4IJiZGeCAYGZmiQOCmZkBDghmZpY4IJiZGeCAYGZmyf8HE+Umv7kCpVAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1134d8910>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# curve fitting for mean of correlation values\n",
    "# Initial guess for curve fitting parameters\n",
    "p0 = [1, 1, 1, 1]\n",
    "# observations\n",
    "x = npoints_list\n",
    "y_meas = avg_corr\n",
    "# least-squares fitting\n",
    "plsq = leastsq(residuals, p0, args=(y_meas, x)) \n",
    "equation = 'y = ((A-D)/(1.0+((x/C)**B))) + D' \n",
    "A = plsq[0][0]\n",
    "B = plsq[0][1]\n",
    "C = plsq[0][2]\n",
    "D = plsq[0][3]\n",
    "\n",
    "# sequence of values for curve plotting\n",
    "xx=np.arange(0,10500,25)\n",
    "# plot fitted curve\n",
    "plt.plot(xx, peval(xx,plsq[0]), color='0.5', lw=0.9)\n",
    "plt.plot(x, y_meas, 'o', color='0.5', ms=5)\n",
    "plt.title(dem)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "clean_rand()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# end GRASS GIS session\n",
    "os.remove(rcfile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
